---
jupytext:
  cell_metadata_json: true
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:init_cell: true

import mmf_setup

mmf_setup.nbinit()
%pylab inline --no-import-all
```

# Travelling Waves: Shooting

+++

Here we consider shooting for solutions to the traveling wave problem.  For more discussion see:

* [Travelling Waves](Travelling%20Waves.ipynb)

+++

## GPE

+++

We start with solutions to the GPE in 1D:

$$
  \I\hbar \dot{\psi}(x, t) =
  -\frac{\hbar^2\psi''(x, t)}{2m}
  + g\abs{\psi(x, t)}^2\psi(x, t).
$$

Implementing a full Galilean transformation, we have

$$
  \psi(x, t) = \sqrt{\rho_V(x_V)}e^{\I[\phi(x_V) - \mu t/\hbar]}e^{\I\bigl[mVx_V + \tfrac{m}{2}V^2t \bigr]/\hbar}
  = \psi_{V}(x_V)e^{\I\bigl[mVx_V + \bigl(\tfrac{m}{2}V^2 - \mu\bigr)t \bigr]/\hbar},\\
  \mu\psi_V(x_V) = -\frac{\hbar^2\psi_V''(x_V)}{2m}
    + g\abs{\psi_V(x_V)}^2\psi_V(x_V),\\
  \psi_V''(x_V) = \frac{2m}{\hbar^2}\left(g\abs{\psi_V(x_V)}^2-\mu\right)\psi_V(x_V).
$$

This is a complex-valued ODE which we can solve by shooing.  Using the global phase invariance and translational invariance, we can assume that at $x_V=0$, the wavefunction is real and has maximum density.  With the parametrization $\psi_V(x) = \sqrt{n(x)}e^{\I x k(x)}$ this means:

$$
  n(x) = n_0 + \tfrac{1}{2}n''(0)x^2 + \order(x^3),\qquad
  k(x) = k_0 + \order(x^2).
$$

*(It might not be obvious that $k'(0) = 0$, but expanding the GPE about $x=0$ and collecting the imaginary parts of the $x^2$ terms in the GPE requires this.)*

$$
  \psi_V(0) = \sqrt{n_0}, \qquad
  \psi_V'(0) = \I k_0\sqrt{n_0},\quad
  \psi_V''(0) = \frac{n_0'' - 2n_0k_0^2}{2\sqrt{n_0}},\\
  \frac{\hbar^2 n_0''}{4mn_0} = \overbrace{\frac{\hbar^2 k_0^2}{2m}}^{E_0} + gn_0 - \mu \leq 0.
$$

Thus, we have a continuous family of solutions with $0 < n_0 < (\mu-E_0)/g$.  As independent parameters we may take $\hbar$, $m$, $g$, specifying units, and $0<n_0$, $0<f=gn_0/(\mu-E_0)<1$, and $k_0$ specifying the solution.

+++

### Numerical Example

+++

Here we numerically solve this equation using the scipy IVP solver, then compare this with our exact solution as implemented in the `gpe.exact_solutions.TravellingWave` class.

```{code-cell}
from scipy.integrate import solve_ivp


class ODE:
    hbar = 1.0
    g = 1.0
    m = 1.0
    periods = 3

    n0 = 1.0
    fraction = 0.9
    k0 = 0.0

    @property
    def mu(self):
        return self.g * self.n0 / self.fraction + self.E0

    @property
    def E0(self):
        return (self.hbar * self.k0) ** 2 / 2 / self.m

    @property
    def y0(self):
        """Return the initial state"""
        psi0 = np.sqrt(self.n0)
        dpsi0 = 1j * self.k0 * psi0
        return (psi0, dpsi0)

    def f(self, x, y):
        psi, dpsi = y
        ddpsi = 2 * self.m * (self.g * np.abs(psi) ** 2 - self.mu) * psi / self.hbar**2
        return (dpsi, ddpsi)

    def event(self, x, y):
        """Termination when dpsi=0 again"""
        if x <= 0:
            return -1
        psi, dpsi = y
        return dpsi.real

    event.terminal = True
    event.direction = -1

    def solve(self, max_step=0.1):
        x_max = 10.0
        ivp = solve_ivp(
            self.f,
            (0, x_max),
            self.y0,
            max_step=max_step,
            # events=[self.event],
        )
        return ivp


ode = ODE()
ode.k0 = 0.3
ode.fraction = 0.9
ode.n0 = 1.0
ode.f(0, ode.y0)
```

```{code-cell}
ivp = ode.solve()
print(ivp.message)
x, psi = ivp.t, ivp.y[0]
plt.figure(figsize=(10, 3))
ax1 = plt.subplot(131)
ax1.plot(x, psi.real)
ax1.plot(x, psi.imag)
ax1.set(xlabel="x", ylabel=r"$\psi$")
ax2 = plt.subplot(132)
ax2.plot(x, abs(psi) ** 2)
ax2.set(xlabel="x", ylabel=r"$n=\psi^2$")
# L = ivp.t_events[0]
```